Aggregating Earth Surface data using GEE and R

Agenda

• From a regular lattice of points create squares
• In each square assign values from aggregating data (e.g. height above sea level of Earth surface)

 

R CRAN

Starters

We use RStudio IDE. Create a .R file (ours is called webinar1_2020_07_29.R) and write the following lines of code. Ctrl-Return to run the line that cursor is on or run all highlighted lines. Or use the button in top right part of the editor.

To add functionalities to your .R script you must load the specific libraries. REMEMBER: make sure that they are installed in your system. If not you can install through RStudio on menu => Tools => Install Packages or directly through R console with command {r eval=F } install.packages("NAME OF PACKAGE")

Below we add three libraries which add functionalities that we need

library(sf)
library(mapview)
library(raster)

 

Coordinate Reference System

My custom Coordinate Reference System (CRS) projection, Lambert Conical Conformal (LCC), NOT secant but tangent at lat=45.827 and lon=11.625 - for more info:

• PROJ library
  
• Wikipedia
  
• Wolfram Mathworld
  
• John P. Snyder “Map projections: A working manual” (pp. 104-110)

myproj <- "+proj=lcc +lat_1=45.827  +lat_2=45.827  +lat_0=45.827 +lon_0=11.625 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

 
 

Read lattice

My points (in lat long) over the regular grid/lattice these points are 666.67 m apart

grid.points <- st_read( "data/webinar1_2020_07_29/points.shp" )

## Reading layer `points' from data source `/archivio/R/shared/webinar/data/webinar1_2020_07_29/points.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1315 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11.83267 ymin: 46.33037 xmax: 12.14322 ymax: 46.55611
## geographic CRS: WGS 84

 

 

Transform CRS

Convert my points from a Geographic (latitude and longitude degrees) to my custom CRS projection

grid.points.myproj <-   grid.points %>% st_transform(myproj)

 

 

View in webmap

Run this line to view the points

  mapview( grid.points.myproj )

 

 

Points to Square

Points to a square Polygon area. How? It is trivial, but faster way was to create a rhomboid with a certain distance from the point (“radius”) using a buffer of half the distance between points (333.33 m) . The minimum bounding box of the buffer with that radius will be a square with sides 2x the size of the radius of the buffer.

NB: nQuadSegs is 1 which means a chord(segment) at each quadrand, i.e. 1/4th of a circunference; this give a rhomboid which is a “very very rough circle” - it saves a lot of memory from drawing with a defaul value of nQuadSegs=30 - it makes a difference when processing a large number of points.

We buffer each feature (point) and view results

nodes.buffered<-st_buffer(grid.points.myproj, 333.33,  nQuadSegs = 1)
mapview( nodes.buffered )

 

 

Then we apply a function to each “rhomboid” for creating the square. It is a bit more complex as the function includes another function inside.

st_bbox_by_feature = function(geom) {
  
  ## make sure that object "geom" becomes a geometry object;
  ## geom is your data set with all the single geometries/polygons
  geom2 = st_geometry(geom)
  
  #' Function to:
  #' (a) take a single geometry, 
  #' (b) create a bounding box with st_bbox and 
  #' (c) convert the bbox object to a new geometry (a square)
  f <- function(single.geom) { 
    st_as_sfc( st_bbox(single.geom,), crs=myproj)
  }
  
  #' This line calls the function above LOOPING over each single geometry (lapply)
  #'  in the geom2 object 
  do.call("c", lapply(geom2, f))
}

This line calls the above function over our data (rhomboids)

tiles <- st_bbox_by_feature( nodes.buffered )

 

Assign the CRS to the new dataset as it got lost. To check the CRS of our data we can use {r} st_crs(tiles)

tiles <- tiles %>% st_set_crs(myproj) 

 

We view all results all over the same map (small subset of data - the first 20 features).

mapview( grid.points.myproj[1:20,] ) + mapview( nodes.buffered[1:20,] )  + mapview( tiles[1:20,] )

 

Save result

Convert to Latitude and Longitude and save result to shapefile for future upload to Google Earth Engine

tiles.latlng <- tiles %>% st_transform("+init=epsg:4326")
#st_write(tiles.latlng, "data/tiles.shp" )

GOOGLE EARTH ENGINE (GEE)

Import shapefile files “tiles.XXX” by left panel “assets” and click “NEW” and choose shapefile from dropdown menu.

DO THINGS in GEE

code shared https://earthengine.googlesource.com/users/2020_Kanan/summer_webinar_series in script https://earthengine.googlesource.com/users/2020_Kanan/summer_webinar_series/+/636ab64efd7b09a0d45d3d10a992b8a44bd119c7/webinar1

EXPORT THE RESULTS TO A SHAPEFILE CALLED tiles_withData.shp

Read the file exported from GEE

grid.points.gee <- st_read( "data/webinar1_2020_07_29/tiles_withData.shp" )

## Reading layer `tiles_withData' from data source `/archivio/R/shared/webinar/data/webinar1_2020_07_29/tiles_withData.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1315 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 11.82833 ymin: 46.32735 xmax: 12.14759 ymax: 46.55912
## geographic CRS: WGS 84

R - Plot

Calculate new IQR field

Create a new attribute column with the interquartile range (The difference between p75 and p25, i.e. 75^th and 25^th percentile)

grid.points.gee$iqr<-grid.points.gee$p75 - grid.points.gee$p25

PLOT variables

pal = mapviewPalette("mapviewSpectralColors")
mapview( grid.points.gee, col.regions = pal(100), zcol = "p50", 
         color=NULL, alpha.regions=0.8  ) 

Figure 1 - 50th Percentile (median)

mapview( grid.points.gee, col.regions = pal(100), zcol = "p10", 
         color=NULL, alpha.regions=0.8  ) 

Figure 2 - 10th Percentile

mapview( grid.points.gee, col.regions = pal(100), zcol = "iqr", 
         color=NULL, alpha.regions=0.8)

Figure 3 - Inter Quartile Range